#include<iostream>
#include<vector>

#include "global.hpp"
#include "predefined.hpp"
#include "solver.hpp"
#include "mesh.hpp"
#include "meanflow.hpp"
#include "mathmatics.hpp"
#ifdef _MPI
#include <mpi.h>
#endif

using namespace std;

void Solver::flux(Mesh &msh, Node *node, Solution *sol, Flux *flux, Meanflow *mf){
    int iZone,iLZone,nNode,iNode,iLNode,iDim,iEq,iVar;
    double Q[NEQ],dQ[NEQ][NDIM],AQx[NEQ],BQy[NEQ],DQ[NEQ];
    double MQ[NEQ],dMQ[NEQ][NDIM];
    #ifndef _2D
    double CQz[NEQ];
    #endif
    vector<vector<double>> A,B,C,D,Vxx,Vxy,Vxz,Vyy,Vyz,Vzz;
    vector<double> x;
    for (iLZone=0;iLZone<msh.nLZone;++iLZone){
        iZone=msh.lZoneList[iLZone];
        nNode=msh.nNode[iZone];
        for (iNode=0;iNode<nNode;++iNode){
            iLNode=msh.idxLNode[iZone]+iNode;
            //msh.calijk(iZone,iNode,ijk);
            for (iEq=0;iEq<NEQ;++iEq){
                for (iDim=0;iDim<NDIM;++iDim){
                    dQ[iEq][iDim]=(flux+iLNode)->dQ[iEq][iDim];
                    dMQ[iEq][iDim]=(mf+iLNode)->dMQ[iEq][iDim];
                }
                Q[iEq]=(sol+iLNode)->Q[iEq];
                MQ[iEq]=(mf+iLNode)->MQ[iEq];
                //if ((dQ[iVar][0]+dQ[iVar][1])*(dQ[iVar][0]+dQ[iVar][1])>TINYNUM){
                
                //}
            }
            //iVar=2;
            //logCtrl.tmpMsg="iVar "+to_string(iVar)+"  Q "+ to_string(Q[iVar])+" dQ "+to_string(dQ[iVar][0]+dQ[iVar][1]);
            //logCtrl.info(logCtrl.tmpMsg.c_str());
            //Calcuate rhs
            m_AdQ(AQx,dQ,MQ);
            m_BdQ(BQy,dQ,MQ);
#           ifndef _2D
            m_CdQ(CQz,dQ,MQ);
#           endif
            m_DQ(DQ,Q,MQ,dMQ);
            #ifndef _2D
            for (iEq=0;iEq<NEQ;++iEq) (flux+iLNode)->FRK[iEq]=-(AQx[iEq]+BQy[iEq]+CQz[iEq]+DQ[iEq]);
            #else
            for (iEq=0;iEq<NEQ;++iEq){
                if ((node+iLNode)->bcTyp>=0) (flux+iLNode)->FRK[iEq]=-(AQx[iEq]+BQy[iEq]+DQ[iEq]);
            } 
            //logCtrl.tmpMsg="iLNode "+to_string(iLNode)+"A "+to_string(A[iEq])+" B "+to_string(B[iEq])+" D "+to_string(D[iEq])+" Q "+to_string((sol+iLNode)->Q[iEq]);
            //logCtrl.info(logCtrl.tmpMsg.c_str());
            #endif
            if ((node+iLNode)->bcTyp<0) {
                for (iEq=0;iEq<NEQ;++iEq) {
                    //(flux+iLNode)->FRK[iEq]=0;
                    //(flux+iLNode)->WRK[iEq]=0;
                }
            }
        }
    }       
}

void Solver::m_AdQ(double *AQx,double dQ[][NDIM],double MQ[NEQ]){
    *(AQx+DENSITY)=MQ[U_X]*dQ[DENSITY][0]+MQ[DENSITY]*dQ[U_X][0];
    *(AQx+U_X)=MQ[U_X]*dQ[U_X][0]+1.0/MQ[DENSITY]*dQ[PRESSURE][0];
    *(AQx+U_Y)=MQ[U_X]*dQ[U_Y][0];
#   ifndef _2D
    *(AQx+U_Z)=MQ[U_X]*dQ[U_Z][0];
#   endif
    *(AQx+PRESSURE)=inputCtrl.gamma*MQ[PRESSURE]*dQ[U_X][0]+MQ[U_X]*dQ[PRESSURE][0];
}

void Solver::m_BdQ(double *BQy,double dQ[][NDIM],double MQ[NEQ]){
    *(BQy+DENSITY)=MQ[U_Y]*dQ[DENSITY][1]+MQ[DENSITY]*dQ[U_Y][1];
    *(BQy+U_X)=MQ[U_Y]*dQ[U_X][1];
    *(BQy+U_Y)=MQ[U_Y]*dQ[U_Y][1]+1.0/MQ[DENSITY]*dQ[PRESSURE][1];
#   ifndef _2D
    *(BQy+U_Z)=MQ[U_Y]*dQ[U_Z][1];
#   endif
    *(BQy+PRESSURE)=inputCtrl.gamma*MQ[PRESSURE]*dQ[U_Y][1]+MQ[U_Y]*dQ[PRESSURE][1];
}


#   ifndef _2D
void Solver::m_CdQ(double *CQz,double dQ[][NDIM],double MQ[NEQ]){
    *(CQz+DENSITY)=MQ[U_Z]*dQ[DENSITY][2]+MQ[0]*dQ[U_Z][2];
    *(CQz+U_X)=MQ[U_Z]*dQ[U_X][2];
    *(CQz+U_Y)=MQ[U_Z]*dQ[U_Y][2];
    *(CQz+U_Z)=MQ[U_Z]*dQ[U_Z][2]+1.0/MQ[0]*dQ[PRESSURE][2];
    *(CQz+PRESSURE)=inputCtrl.gamma*MQ[PRESSURE]*dQ[U_Z][2]+MQ[U_Z]*dQ[PRESSURE][2];
}
#   endif

void Solver::m_DQ(double *DQ,double Q[NEQ],double MQ[NEQ],double dMQ[][NDIM]){
    int iDim,iEq,i;
    double divMU=0,UgradU[NEQ],QgradMQ[NEQ];
    for (iDim=0;iDim<NDIM;++iDim){
        divMU+=dMQ[iDim+1][iDim];
        UgradU[iDim+1]=0;
        for (i=0;i<NDIM;++i)UgradU[iDim+1]+=MQ[iDim+1]*dMQ[iDim+1][i];
        UgradU[iDim+1]=UgradU[iDim+1]/MQ[DENSITY]*Q[DENSITY];
    }
    UgradU[DENSITY]=divMU*Q[DENSITY];
    UgradU[PRESSURE]=0;
    for (iEq=0;iEq<NEQ;++iEq){
        QgradMQ[iEq]=0;
        for (iDim=0;iDim<NDIM;++iDim) QgradMQ[iEq]+=Q[iDim+1]*dMQ[iEq][iDim];
    }
    for (iEq=0;iEq<NEQ;++iEq) *(DQ+iEq)=UgradU[iEq]+QgradMQ[iEq];
    *(DQ+PRESSURE)+=inputCtrl.gamma*Q[PRESSURE]*divMU;
}

void Solver::RK(int iRK,int seRK,Mesh &msh, Solution *sol, Flux *flux){
    int iLNode,iEq;
    double dt=inputCtrl.dt;
    for (iLNode=0;iLNode<msh.nLNode;++iLNode){
        for (iEq=0;iEq<NEQ;++iEq){
            (flux+iLNode)->WRK[iEq]=s.rk[seRK][0][iRK]*((flux+iLNode)->WRK[iEq])+dt*((flux+iLNode)->FRK[iEq]);
            (sol+iLNode)->Q[iEq]+=s.rk[seRK][1][iRK]*((flux+iLNode)->WRK[iEq]);
            //logCtrl.tmpMsg="iLNode "+to_string(iLNode)+"WRK "+to_string((flux+iLNode)->WRK[iEq])+" FRK "+to_string((flux+iLNode)->FRK[iEq])+" Q "+to_string((sol+iLNode)->Q[iEq]);
            //logCtrl.info(logCtrl.tmpMsg.c_str());
        }
    }
}